Series de tiempo¶
De entre las series de tiempo para cada enfermedad en el periodo descrito, no todas contienen información de casos reportados durante el periodo de tiempo establecido para la investigación, por lo que se extraen los que reportan al menos la mitad del periodo (260 semanas)
In [142]:
# https://stackoverflow.com/a/16916611
print('Iniciales {}'.format(len(cie)))
cie = cie.filter(lambda x: x['sem'].count() >= 260 )
cie.reset_index(drop=True, inplace=True)
cie = cie.groupby('cie')
print('Restantes {}'.format(len(cie)))
In [143]:
print('Iniciales {}'.format(len(cieG)))
cieG = cieG.filter(lambda x: x['sem'].count() >= 260 )
cieG.reset_index(drop=True, inplace=True)
cieG = cieG.groupby(cieG.cie.str[0])
print('Restantes {}'.format(len(cieG)))
AsÃ, de 138 series de tiempo de enfermedades, se obtienen 40 en las que al menos se cuenta con datos semanales de 5 años. Para dichas enfermedades se obtienen los pesos de la regresión lineal y se obtiene la serie de tiempo sin la tendencia y las autocorrelaciones (eliminando la a92.3 porque viene vacÃa)
In [490]:
from scipy import signal
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf
from statsmodels.tsa.stattools import acf
ciesF = [] # CIEs CaracterÃsticas
ciesTSt = [] # CIEs Series de tiempo
for name, group in cie:
if name == 'a92.3':
continue
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html
detrended = signal.detrend(group.casos)
a, b, r, p, e = stats.linregress(group['sem'], group.casos)
print("y = f(x) = {} x + {}".format(a, b))
print("error", e)
print("p = ", p)
print("pendiente {:s}significativa".format("no " if p >= 0.05 else ""))
print("R^2", r**2)
plt.figure(figsize=(12, 2))
plt.plot(group['sem'], group.casos)
plt.plot(group['sem'], detrended, c='black')
plt.plot(group['sem'], (a * group['sem'] + b), label = 'y = {:.1f}x + {:.0f}'.format(a, b), color = 'red', linewidth = 3)
plt.title(name)
plt.xlabel("Semana")
plt.ylabel("Casos normalizados")
plt.show()
# https://stackoverflow.com/questions/48497756/time-series-distance-metric
plt.figure(figsize=(12, 2))
# https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.cumsum.html
plt.plot(group['sem'], group.casos.cumsum(), c='green')
plt.title(name)
plt.xlabel("Semana")
plt.ylabel("Acumulado de Casos normalizados")
plt.show()
# https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
plot_acf(detrended, lags=52)
# https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.acf.html
plt.title(name)
plt.xlabel("Retraso en semanas")
plt.ylabel('Correlación')
plt.show()
temp = [a, b]
# https://stackoverflow.com/a/3748071
temp.extend(acf(detrended, nlags=52))
temp.append(name)
ciesF.append(temp)
temp2 = list(group.casos)
#temp2.append(name)
ciesTSt.append(temp2)
In [486]:
for name, group in cieG:
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html
detrended = signal.detrend(group.casos)
a, b, r, p, e = stats.linregress(group['sem'], group.casos)
print("y = f(x) = {} x + {}".format(a, b))
print("error", e)
print("p = ", p)
print("pendiente {:s}significativa".format("no " if p >= 0.05 else ""))
print("R^2", r**2)
plt.figure(figsize=(12, 2))
plt.plot(group['sem'], group.casos)
plt.plot(group['sem'], detrended, c='black')
plt.plot(group['sem'], (a * group['sem'] + b), label = 'y = {:.1f}x + {:.0f}'.format(a, b), color = 'red', linewidth = 3)
plt.title(name)
plt.xlabel("Semana")
plt.ylabel("Casos normalizados")
plt.show()
# https://stackoverflow.com/questions/48497756/time-series-distance-metric
plt.figure(figsize=(12, 2))
# https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.cumsum.html
plt.plot(group['sem'], group.casos.cumsum(), c='green')
plt.title(name)
plt.xlabel("Semana")
plt.ylabel("Acumulado de Casos normalizados")
plt.show()
plt.figure()
# https://machinelearningmastery.com/gentle-introduction-autocorrelation-partial-autocorrelation/
plot_acf(detrended, lags=52)
# https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.acf.html
plt.title(name)
plt.xlabel("Retraso en semanas")
plt.ylabel('Correlación')
plt.show()
Se extraen las caracterÃsticas de cada CIE en tanto serie de tiempo. A saber, su pendiente, ordenada en el origen y las autocorrelaciones con retraso de 1 a 52 semanas (eliminando el retraso de 0 semanas)
In [491]:
ciesF = pd.DataFrame(ciesF)
# https://stackoverflow.com/a/11346337
colNames = ['m', 'b']
for i in range (53):
colNames.append('ac' + str(i))
colNames.append('cie')
ciesF.columns = colNames
ciesF = ciesF.drop(['ac0'], axis=1)
ciesF.sort_values(by=['m'], ascending=False)
Out[491]:
Se agrupan las CIEs por casos dados en una semana, como series de tiempo.
In [492]:
ciesTS = pd.DataFrame(ciesTSt)
t = list(cie.groups.keys())
t.remove('a92.3')
ciesTS['cie'] = t
ciesTS.sample(3)
Out[492]:
Pronóstico¶
Una vez eliminada la tendencia se puede comprobar que las series de tiempo para cada enfermedad son estacionales
In [493]:
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries, w, name):
#Determing rolling statistics
rolmean = timeseries.rolling(w).mean()
rolstd = timeseries.rolling(w).std()
#Plot rolling statistics:
orig = plt.plot(timeseries,label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean')
std = plt.plot(rolstd, color='black', label = 'Rolling Std')
plt.legend(loc='best')
plt.title(name)
plt.show(block=False)
#Perform Dickey-Fuller test:
print('Results of Dickey-Fuller Test:')
dftest = adfuller(timeseries, autolag='AIC')
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
print(dfoutput)
for i in range(len(ciesTS)):
temp = ciesTS.iloc[i, : 508]
temp.dropna(inplace=True)
test_stationarity(temp, 52, ciesTS.iloc[i, -1])
Por ello es posible pronosticarlas con el método de Holt-Winter
In [494]:
# https://www.analyticsvidhya.com/blog/2018/02/time-series-forecasting-methods/
# https://machinelearningmastery.com/time-series-forecasting-methods-in-python-cheat-sheet/
# https://towardsdatascience.com/time-series-in-python-exponential-smoothing-and-arima-processes-2c67f2a52788
from statsmodels.tsa.api import Holt, SimpleExpSmoothing, ExponentialSmoothing
pronosticos = []
real = []
for i in range(len(ciesTS)):
temp = ciesTS.iloc[i, : 508]
temp.dropna(inplace=True)
# https://medium.com/datadriveninvestor/how-to-build-exponential-smoothing-models-using-python-simple-exponential-smoothing-holt-and-da371189e1a1
# Train = 0.7
train = round(0.7 * len(temp))
f = ExponentialSmoothing(np.asarray(temp.iloc[0:train])).fit(smoothing_level = 0.1)
# https://stackoverflow.com/a/50786171
pred = f.predict(start=train + 1, end=len(temp))
fcast = f.forecast(len(temp) - train)
plt.figure(figsize=(12, 4))
plt.plot(temp)
plt.plot(f.fittedvalues, c='black')
plt.plot(range(train, len(temp)), fcast, c='red')
#plt.plot(range(train, len(temp)), pred, c='red')
plt.title(ciesTS.iloc[i, 508])
plt.legend(["Serie de tiempo", "Ajuste de Holt-Winter", "Pronóstico de Holt-Winter"])
plt.show()
# https://stackoverflow.com/a/15863028
real.append(temp.iloc[-1])
pronosticos.append(fcast[0])
Que presenta muy buen ajuste respecto a los datos reales
In [417]:
a, b, r, p, e = stats.linregress(real, pronosticos)
print("y = f(x) = {:.4f} x + {:.4f}".format(a, b))
print("error", e)
print("p = ", p)
print("pendiente {:s}significativa".format("no " if p >= 0.05 else ""))
print("R^2", r**2)
plt.title('Logaritmo del pronóstico', fontsize = 20)
plt.xlabel('Logaritmo del último caso normalizado')
plt.ylabel('Pronóstico')
plt.scatter(real, pronosticos)
plt.show()
Lo que se evidencia claramente al utilizar escala logarÃtmica en los datos
In [536]:
plt.title('Logaritmo del pronóstico', fontsize = 20)
plt.xlabel('Logaritmo del último caso normalizado')
plt.ylabel('Pronóstico')
plt.scatter(np.log(real), np.log(pronosticos))
plt.show()
In [526]:
# https://machinelearningmastery.com/time-series-forecast-uncertainty-using-confidence-intervals-python/
# https://machinelearningmastery.com/make-sample-forecasts-arima-python/
from statsmodels.tsa.arima_model import ARIMA
def difference(dataset, interval=1):
diff = list()
for i in range(interval, len(dataset)):
value = dataset[i] - dataset[i - interval]
diff.append(value)
return pd.Series(diff)
for i in range(len(ciesTS)):
temp = ciesTS.iloc[i, : 508]
temp.dropna(inplace=True)
train = round(0.9 * len(temp))
ts = difference(np.asarray(temp), interval=len(temp) - train - 1)
m = ARIMA(ts, order=(3, 1, 0))
f = m.fit(trend = 'nc', disp = 0)
pred = f.predict(start=train + 1, end=len(temp))
tempDiff = temp - temp.shift()
plt.figure(figsize=(16, 4))
plt.plot(tempDiff)
plt.plot(f.fittedvalues, color='red')
plt.plot(range(train, len(temp)), pred, c='black')
plt.legend(['Diferencia', 'Valores ajustados de ARIMA', 'Pronóstico de ARIMA'])
plt.show()
